#####################################################################
# Public Preferences for Intergroup Assistance in Conflicts Facing  #
# Joint External Threats: Lessons from COVID-19 in Israel           #
# Liran Harsgor and Alon Yakter                                     #
# Journal of Conflict Resolition                                    #
#                                                                   #
# Replication code - Conjoint analysis                              #
#####################################################################
# To load the data, users should place the relevant csv	file in     #
# in their R working directory (or type the full file paths as      #
# saved on their computer).                                         #
#####################################################################

## Required packages
library(cjoint)
library(foreign)
library(dplyr)
library(cregg)
library(ggplot2)
library(patchwork)
library(labelled)
library(stargazer)

## set-up data 
# load data file
dat.bi <- read.csv("Conjoint experiment july.csv")
# clean empty responses
dat.bi.na <- dat.bi %>% filter(!is.na(selected))
# set trait category order for the plot
dat.bi.na$action <- factor(dat.bi.na$action, levels = c("Lockdown","Worker Ban","Monitor","Protective Equipment","Medical Aid"))
dat.bi.na$effect_pal <- factor(dat.bi.na$effect_pal, levels = c("Deterioration","No Effect","Improvement"))
dat.bi.na$effect_isr <- factor(dat.bi.na$effect_isr, levels = c("Fewer Infections","No Change"))
dat.bi.na$funding <- factor(dat.bi.na$funding, levels = c("Palestinian Taxes","Half-and-Half","Israeli Budget"))
dat.bi.na$collaboration <- factor(dat.bi.na$collaboration, levels = c("No Coordination","Only With PA","With PA and Hamas"))
# set variable labels
var_label(dat.bi.na) = list(action = "POLICY TYPE", effect_pal = "PALESTINIAN ILLNESS", effect_isr = "CROSS-INFECTIONS",
               funding = "FUNDING SOURCE", collaboration = "COORDINATION")

##########################
# full sample (figure 2) #
##########################
# analysis
out.all <- cj(data = dat.bi.na,
          selected ~ action * effect_pal + effect_isr + funding + collaboration ,
          id = ~Response.ID,
          estimate = "mm",
)
# plot outcomes
cj.all <- plot(out.all, vline = 0.5) + 
  ggplot2::theme(legend.position = "none") 
ggsave(file = "Figure 2.pdf", plot=cj.all, width=4, height=5)

############################################################
# heterogeneity by health and economic concerns (Figure 3) #
############################################################
## health concerns
# split sample by low, medium, high concerns
dat.bi.na$health_concern <- as.numeric(dat.bi.na$health_concern)
dat.bi.na$conc.h <- cut(dat.bi.na$health_concern,c(0,2,3,5))
# analysis
out.h <- cj(data = dat.bi.na,
            selected ~ action * effect_pal + effect_isr + funding + collaboration ,
            id = ~Response.ID,
            estimate = "mm",
            by=~conc.h,
)
# F-test for intergroup differences
cj_anova(data = na.omit(dat.bi.na,na.action="omit"),
         selected ~ action * effect_pal + effect_isr + funding + collaboration ,
         by=~conc.h,
)
# plot outcomes (left-side panel)
cj.health <- plot(out.h, group = "conc.h", vline = 0.5, legend_title = "Health Concerns") +
  ggplot2::scale_color_manual(labels=c("Low","Medium","High",""),values=c("#1B1919FF","#f37735","#00b159")) +
  ggplot2::theme(legend.title=element_blank()) + 
  ggplot2::ggtitle("Health Concern")

## economic concerns
# split sample by low, medium, high concerns
dat.bi.na$econ_concern <- as.numeric(dat.bi.na$econ_concern)
dat.bi.na$conc.e <- cut(dat.bi.na$econ_concern,c(0,2,3,5))
# analysis
out.e <- cj(data = dat.bi.na,
            selected ~ action * effect_pal + effect_isr + funding + collaboration ,
            id = ~Response.ID,
            estimate = "mm",
            by=~conc.e,
)
# F-test for intergroup differences
cj_anova(data = na.omit(dat.bi.na,na.action="omit"),
         selected ~ action * effect_pal + effect_isr + funding + collaboration ,
         by=~conc.e,
)
# plot outcomes (right-hand panel)
cj.econ <- plot(out.e, group = "conc.e", vline = 0.5, legend_title = "Economic Concerns") +
  ggplot2::scale_color_manual(labels=c("Low","Medium","High",""),values=c("#1B1919FF","#f37735","#42B540FF")) +
  ggplot2::theme(legend.title=element_blank()) + 
  ggplot2::ggtitle("Economic Concern")

## combined graph (figure 3)
plot.h.e <- cj.health + cj.econ
ggsave(file = "Figure 3.pdf", plot=plot.h.e, width=7, height=5.5)

#########################################################
# heterogeneity by left-right identification (figure 4) #
#########################################################
# split sample by left, center, right
dat.bi.na$par <- cut(dat.bi.na$left_right,c(0,3,4,7))
dat.bi.na$par <- factor(dat.bi.na$par, levels = c("(4,7]", "(3,4]", "(0,3]"))
# analysis
out.lr <- cj(data = dat.bi.na,
          selected ~ action * effect_pal + effect_isr + funding + collaboration ,
          id = ~Response.ID,
          estimate = "mm",
          by=~par,
)
# F-test for intergroup differences
cj_anova(data = na.omit(dat.bi.na,na.action="omit"),
         selected ~ action * effect_pal + effect_isr + funding + collaboration ,
         by=~par,
)
# plot outcomes
cj.lr <- plot(out.lr, group = "par", vline = 0.5, legend_title = "Ideology") +
  ggplot2::scale_color_manual(labels=c("Right","Center","Left",""),values=c("#619CFF","#7CAE00","#F8766D")) + 
  ggplot2::theme(legend.title=element_blank()) 
ggsave(file = "Figure 3.pdf", plot=cj.lr, width=4, height=5.5)

################################################
# presentation of full results (tables A4-A11) #
################################################
## full sample
#  marginal means (table A4)
t.all <- head(out.all[c("feature", "level", "estimate", "std.error")],60L)
stargazer(t.all, summary=FALSE, rownames = FALSE)
# full sample - ACME (table A5)
out.all.amce <- cj(data = dat.bi.na,
                   selected ~ action * effect_pal + effect_isr + funding + collaboration ,
                   id = ~Response.ID,
)
t.all.amce <- head(out.all.amce[c("feature", "level", "estimate", "std.error")],60L)
stargazer(t.all.amce, summary=FALSE, rownames = FALSE)
## subgroup analysis by ideology 
# marginal means (table A6) 
out.lr <- head(out.lr[c("feature", "level", "estimate", "std.error")],60L)
stargazer(t.lr, summary=FALSE, rownames = FALSE)
# ACME (table A7) 
out.lr.amce <- cj(data = dat.bi.na,
             selected ~ action * effect_pal + effect_isr + funding + collaboration ,
             id = ~Response.ID,
             by=~par,
)
out.lr.amce <- head(out.lr.amce[c("feature", "level", "estimate", "std.error")],60L)
stargazer(out.lr.amce, summary=FALSE, rownames = FALSE)

## subgroup analysis by health concerns 
# marginal means (table A8) 
t.h <- head(out.h[c("feature", "level", "estimate", "std.error")],60L)
stargazer(t.h, summary=FALSE, rownames = FALSE)
# ACME (table A9) 
out.h.amce <- cj(data = dat.bi.na,
                 selected ~ action * effect_pal + effect_isr + funding + collaboration ,
                 id = ~Response.ID,
                 by=~conc.h,
)
t.h.amce <- head(out.h.amce[c("feature", "level", "estimate", "std.error")],60L)
stargazer(t.h.amce, summary=FALSE, rownames = FALSE)

## subgroup analysis by economic concerns 
# marginal means (table A10) 
t.e <- head(out.e[c("feature", "level", "estimate", "std.error")],60L)
stargazer(t.e, summary=FALSE, rownames = FALSE)
# ACME (table A11) 
out.e.amce <- cj(data = dat.bi.na,
                 selected ~ action * effect_pal + effect_isr + funding + collaboration ,
                 id = ~Response.ID,
                 by=~conc.e,
)
t.e.amce <- head(out.e.amce[c("feature", "level", "estimate", "std.error")],60L)
stargazer(t.e.amce, summary=FALSE, rownames = FALSE)

############################################
# robustness: interactions (Figures A3-A4) #
############################################
### heterogeneity by health concerns and left-right (figure A3)
## heterogeneity by health concerns among lefitsts
# analysis
out.lr.h.l <- cj(data = filter(dat.bi.na,left_right<=3),
             selected ~ action * effect_pal + effect_isr + funding + collaboration ,
             id = ~Response.ID,
             estimate = "mm",
             by=~conc.h,
)
# plot 
cj.lr.h.l <- plot(out.lr.h.l, group = "conc.h", vline = 0.5, legend_title = "Health Concerns") +
  ggplot2::scale_color_manual(labels=c("Low","Medium","High",""),values=c("#1B1919FF","#f37735","#42B540FF")) +
  ggplot2::theme(legend.title=element_blank()) +
  ggplot2::ggtitle("Leftists")
## heterogeneity by health concerns among centrists
# analysis
out.lr.h.c <- cj(data = filter(dat.bi.na,left_right==4),
                 selected ~ action * effect_pal + effect_isr + funding + collaboration ,
                 id = ~Response.ID,
                 estimate = "mm",
                 by=~conc.h,
)
# plot
cj.lr.h.c <- plot(out.lr.h.c, group = "conc.h", vline = 0.5, legend_title = "Health Concerns") +
  ggplot2::scale_color_manual(labels=c("Low","Medium","High",""),values=c("#1B1919FF","#f37735","#42B540FF")) +
  ggplot2::theme(legend.title=element_blank()) +
  ggplot2::ggtitle("Centrists")
## heterogeneity by health concerns among rightists
# analysis
out.lr.h.r <- cj(data = filter(dat.bi.na,left_right>=5),
                 selected ~ action * effect_pal + effect_isr + funding + collaboration ,
                 id = ~Response.ID,
                 estimate = "mm",
                 by=~conc.h,
)
# plot
cj.lr.h.r <- plot(out.lr.h.r, group = "conc.h", vline = 0.5, legend_title = "Health Concerns") +
  ggplot2::scale_color_manual(labels=c("Low","Medium","High",""),values=c("#1B1919FF","#f37735","#42B540FF")) +
  ggplot2::theme(legend.title=element_blank()) +
  ggplot2::ggtitle("Rightists") 
## combine graphs
plot.lr.h <- cj.lr.h.l + cj.lr.h.c + cj.lr.h.r + 
  plot_annotation(title = "Health Concerns by Left-Right Identification")
ggsave(file = "Figure A3.pdf", plot=plot.lr.h, width=9, height=5.5)

### interaction: heterogeneity by  economic concerns and left-right (figure A4)
## heterogeneity by economic concerns among lefitsts
# analysis
out.lr.e.l <- cj(data = filter(dat.bi.na,left_right<=3),
                 selected ~ action * effect_pal + effect_isr + funding + collaboration ,
                 id = ~Response.ID,
                 estimate = "mm",
                 by=~conc.e,
)
# plot 
cj.lr.e.l <- plot(out.lr.e.l, group = "conc.e", vline = 0.5, legend_title = "Health Concerns") +
  ggplot2::scale_color_manual(labels=c("Low","Medium","High",""),values=c("#1B1919FF","#f37735","#42B540FF")) +
  ggplot2::theme(legend.title=element_blank()) +
  ggplot2::ggtitle("Leftists")
## heterogeneity by economic concerns among centrists
# anaysis
out.lr.e.c <- cj(data = filter(dat.bi.na,left_right==4),
                 selected ~ action * effect_pal + effect_isr + funding + collaboration ,
                 id = ~Response.ID,
                 estimate = "mm",
                 by=~conc.e,
)
# plot
cj.lr.e.c <- plot(out.lr.e.c, group = "conc.e", vline = 0.5, legend_title = "Health Concerns") +
  ggplot2::scale_color_manual(labels=c("Low","Medium","High",""),values=c("#1B1919FF","#f37735","#42B540FF")) +
  ggplot2::theme(legend.title=element_blank()) +
  ggplot2::ggtitle("Centrists")
## heterogeneity by economic concerns among rightists
# analysis
out.lr.e.r <- cj(data = filter(dat.bi.na,left_right>=5),
                 selected ~ action * effect_pal + effect_isr + funding + collaboration ,
                 id = ~Response.ID,
                 estimate = "mm",
                 by=~conc.e,
)
# plot
cj.lr.e.r <- plot(out.lr.e.r, group = "conc.e", vline = 0.5, legend_title = "Health Concerns") +
  ggplot2::scale_color_manual(labels=c("Low","Medium","High",""),values=c("#1B1919FF","#f37735","#42B540FF")) +
  ggplot2::theme(legend.title=element_blank()) +
  ggplot2::ggtitle("Rightists") 
## combine graphs
plot.lr.e <- cj.lr.e.l + cj.lr.e.c + cj.lr.e.r + 
  plot_annotation(title = "Health Concerns by Left-Right Identification")
ggsave(file = "Figure A4.pdf", plot=plot.lr.e, width=9, height=5.5)

########################################
# Conjoint diagnostics (figures A5-A8) #
########################################
## Display Frequencies and Proportions (figure A5)
dgnst.freq <- plot(cj_freqs(data = dat.bi.na,
              selected ~ action * effect_pal + effect_isr + funding + collaboration,
              id = ~Response.ID)) + 
  ggplot2::theme(legend.position = "none")
ggsave(file = "Figure A5.pdf", plot=dgnst.freq, width=6, height=5)
# verifying which profiles were constrained from showing (as planned)
subset(cj_props(data = dat.bi.na,
                ~ action + effect_pal, id = ~Response.ID), Proportion == 0)

## Carryover between choice tasks and profiles (figure A6)
# choice tasks (left-side panel)
dat.bi.na$task_factor <- factor(dat.bi.na$task)
dgnst.carry.ts <- cj(data = dat.bi.na,
              selected ~ action * effect_pal + effect_isr + funding + collaboration ,
              id = ~Response.ID,
              estimate = "amce",
              by = ~task_factor,)
dgnst.carry.plot.ts <- plot(dgnst.carry.ts, vline = 0, group = "task_factor") + 
  ggplot2::guides(color=guide_legend(title="Task Order")) + 
  ggplot2::ggtitle("Task Order")
# profiles (right-side panel)
dat.bi.na$profile_factor <- factor(dat.bi.na$profile)
dgnst.carry.pr <- cj(data = dat.bi.na,
                  selected ~ action * effect_pal + effect_isr + funding + collaboration ,
                  id = ~Response.ID,
                  estimate = "amce",
                  by = ~profile_factor,)
dgnst.carry.plot.pr <- plot(dgnst.carry.pr, vline = 0, group = "profile_factor") + 
  ggplot2::guides(color=guide_legend(title="Profile Order")) + 
  ggplot2::ggtitle("Profile Order")
# combined carryover plot
dgnst.carry.plot <- dgnst.carry.plot.ts + dgnst.carry.plot.pr
ggsave(file = "Figure A6.pdf", plot=dgnst.carry.plot, width=7, height=5.5)

## balance testing - demographics (figure A7)
blnce.sex <- plot(mm(data = dat.bi.na,
                     sex ~ action * effect_pal + effect_isr + funding + collaboration,
                     id = ~Response.ID), 
                  vline = mean(dat.bi.na$sex, na.rm = TRUE)) +
  ggplot2::ggtitle("Sex") + 
  ggplot2::theme(legend.position = "none") + 
  ggplot2::theme(axis.text.x = element_blank())
blnce.age <- plot(mm(data = dat.bi.na,
                     age ~ action * effect_pal + effect_isr + funding + collaboration,
                     id = ~Response.ID), 
                  vline = mean(dat.bi.na$age, na.rm = TRUE)) +
  ggplot2::ggtitle("Age") + 
  ggplot2::theme(legend.position = "none") + 
  ggplot2::theme(axis.text = element_blank())
blnce.edu <- plot(mm(data = dat.bi.na,
                     education ~ action * effect_pal + effect_isr + funding + collaboration,
                     id = ~Response.ID), 
                  vline = mean(dat.bi.na$education, na.rm = TRUE)) +
  ggplot2::ggtitle("Education") + 
  ggplot2::theme(legend.position = "none") + 
  ggplot2::theme(axis.text = element_blank())
blnce.relig <- plot(mm(data = dat.bi.na,
                       religiosity ~ action * effect_pal + effect_isr + funding + collaboration,
                       id = ~Response.ID), 
                    vline = mean(dat.bi.na$religiosity, na.rm = TRUE)) +
  ggplot2::ggtitle("Religiosity") + 
  ggplot2::theme(legend.position = "none") + 
  ggplot2::theme(axis.text = element_blank())
blnce.income <- plot(mm(data = dat.bi.na,
                        income ~ action * effect_pal + effect_isr + funding + collaboration,
                        id = ~Response.ID), 
                     vline = mean(dat.bi.na$income, na.rm = TRUE)) +
  ggplot2::ggtitle("Income") + 
  ggplot2::theme(legend.position = "none") + 
  ggplot2::theme(axis.text = element_blank())
blnce.region <- plot(mm(data = dat.bi.na,
                        region ~ action * effect_pal + effect_isr + funding + collaboration,
                        id = ~Response.ID), 
                     vline = mean(dat.bi.na$region, na.rm = TRUE)) +
  ggplot2::ggtitle("Region") + 
  ggplot2::theme(legend.position = "none") + 
  ggplot2::theme(axis.text = element_blank())
blnce.plot.dem <- blnce.sex + blnce.age + blnce.edu + blnce.relig + blnce.income + blnce.region + 
  plot_layout(nrow = 1)
ggsave(file = "Figure A7.pdf", plot=blnce.plot.dem, width=9, height=5.5)

## balance testing - health and economic concerns and ideology (figure A8)
blnce.lr <- plot(mm(data = dat.bi.na,
        left_right ~ action * effect_pal + effect_isr + funding + collaboration,
        id = ~Response.ID), 
     vline = mean(dat.bi.na$left_right, na.rm = TRUE)) +
  ggplot2::ggtitle("Left-Right") + 
  ggplot2::theme(legend.position = "none") + 
  ggplot2::theme(axis.text.x = element_blank())
blnce.h <- plot(mm(data = dat.bi.na,
        health_concern ~ action * effect_pal + effect_isr + funding + collaboration,
        id = ~Response.ID), 
     vline = mean(dat.bi.na$health_concern, na.rm = TRUE)) +
  ggplot2::ggtitle("Health Concern") + 
  ggplot2::theme(legend.position = "none") + 
  ggplot2::theme(axis.text = element_blank())
blnce.e <- plot(mm(data = dat.bi.na,
        econ_concern ~ action * effect_pal + effect_isr + funding + collaboration,
        id = ~Response.ID), 
     vline = mean(dat.bi.na$econ_concern, na.rm = TRUE)) +
  ggplot2::ggtitle("Economic Concern") + 
  ggplot2::theme(legend.position = "none") + 
  ggplot2::theme(axis.text = element_blank())
blnce.plot.iv <- blnce.lr + blnce.h + blnce.e  + 
  plot_layout(nrow = 1)
ggsave(file = "Figure A8.pdf", plot=blnce.plot.iv, width=8, height=5.5)